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Abstract 



I derive precise results for absorbing-state phase transitions using exact (numerically determined) 

quasistationary probability distributions for small systems. Analysis of the contact process on rings 

of 23 or fewer sites yields critical properties (control parameter, order-parameter ratios, and critical 
•4— > . 

exponents z and (3/v±) with an accuracy of better than 0.1%; for the exponent v± the accuracy is 

03 



about 0.5%. Good results are also obtained for the pair contact process. 
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INTRODUCTION 



Stochastic processes with an absorbing state arise frequently in statistical physics 



HQ. 



In systems with spatial structure, phase transitions to an absorbing state, as exemplified 
by the contact process 3|, |4(, are widely studied in connection with self-organized criticality 
5|, the transition to turbulence 6|, and issues of universality in nonequilibrium critical 
phenomena [?], S, S, III]]. Interest in such transitions should continue to grow in the wake 
of experimental confirmation in a liquid crystal system This Letter presents a new 

theoretical approach to absorbing-state phase transitions via analysis of exact (numerical) 
quasistationary (QS) probability distributions. 

The quasistationary probability distribution (QSD) provides a wealth of information 
about systems exhibiting an absorbing-state phase transition [3, [3]. (Since the only true 
stationary state for a finite system is the absorbing one, "stationary-state" simulations in 
fact probe QS properties, that is, conditioned on survival.) In particular, the order parame- 
ter and its moments, static correlation functions, and the QS lifetime are all accessible from 
the QSD. Until now, QS properties of systems with spatial structure have been determined 



only via simulation 13|, uA, ll5j] ; here I develop an effective scheme for determining the QSD 
on rings of L sites. 

The QSD is defined as follows. Consider a continuous-time Markov process X t with 
state A absorbing: if X t = A, then X t > = A at all subsequent times. The transition rates 
Wc,c (from state C to state C) are such that wc,a = 0, VC. (Some processes have several 
absorbing states, A\, A n .) Let pc{t) denote the probability of state C at time t, given 
some initial state Xq ^ A. The survival probability P s {t) = Ylc^APc{t) is the probability 
that the process has not visited the absorbing state up to time t. We suppose that as t — > oo 
the pc, normalized by the survival probability, attain a time-independent form, allowing us 
to define the QSD: 

with p A = 0; it is normalized so: J2c^aPc = 1- 

In principle, one could integrate the master equation numerically and extract the QSD 
in the long-time limit. Such an approach is very time-consuming, and essentially useless 
for processes with a large state space. I use instead the iterative scheme demonstrated in 



[161 ] . Given some initial guess for the distribution p c , the following relation is iterated until 
convergence is achieved: 



V'c = ap c + (1 - a) TC . (2) 
w c - r A 

Here r c = Ylc w c,C'Pc * s the probability flux (in the master equation) into state C {ta 
is the flux to the absorbing state, so that l/r^ gives the lifetime of the QS state), and 
wc = Ylc> w c,c is the total rate of transitions out of state C. The parameter a can take 
any value between and 1; in practice rapid convergence is obtained with a = 0.1. 

The iterative scheme is used to determine the QSD of the contact process (CP) on rings 
of L sites. In the CP , each site i of a lattice is either occupied {o~i(t) = 1), or vacant 

(o~i(t) = 0). Transitions from <ji = 1 to cr, = occur at a rate of unity, independent of the 
neighboring sites. The reverse transition is only possible if at least one neighbor is occupied: 
the transition from o~i = to o~i = 1 occurs at rate Ar, where r is the fraction of nearest 
neighbors of site % that are occupied; thus the state o~i = for all i is absorbing. (A is a 
control parameter governing the rate of spread of activity.) 

Although no exact results are available, the CP has been studied intensively via series 
expansion and Monte Carlo simulation. Since its scaling properties have been discussed 
, |9| we review them only briefly here. The best estimate for the critical point 



extensively [7|, _ 

in one dimension is A c = 3.297848(20), as determined via series analysis [l7]. Approaching 
the critical point, the correlation length £ and lifetime r diverge, following £ cx |A| _I/± and 
r oc |A| _i/ ii, where A = (A — A c )/A c is the relative distance from the critical point. The 
order parameter (the fraction of active sites), scales as p oc A' 3 for A > 0. Near the critical 
point, finite-size scaling (FSS) [181 ]. implies that average properties such as p depend on L 
through the scaling variable AL 1//iy± , leading, at the critical point, to r oc L z , with dynamic 
exponent z = v\\/v±, and p oc L~P' U± . 

The computational algorithm for determining the QSD consists of three components. 
The first (applicable to any model having two states per site) enumerates all configurations 
on a ring of L sites. Configurations differing only by a lattice translation are treated as 
equivalent. In subsequent stages only one representative of each equivalence class is used, 
yielding a considerable speedup and reduction in memory requirements. (The multiplicity, 
or number of configurations associated with each equivalence class, is needed for calculating 
observables.) The second component runs through the list of configurations, enumerating 



all possible transitions. Here proper care must be taken to determine the weight of each 
transition, due to the varying multiplicity of initial and final configurations. The exit rate 
for each configuration C is simply: wc = ric + (A/2) Cc, where n c is the number of occupied 
sites and cq the number of occupied-vacant nearest-neighbor pairs. To determine tq one 
enumerates, for each configuration, all transitions from some other state C to C. (Each 
vacant site i m C implies a transition from a configuration C, differing from C only in 
that site % is occupied; each nearest-neighbor pair of occupied sites + l in C implies 
transitions from a C in which either % or % + 1 is vacant. Transitions between the same pair 
of configurations C and C are grouped together, with the proper multiplicity stored in the 
associated weight.) The final part of the algorithm determines the QSD via the iterative 
procedure described above. The specific rules of the model enter only in the second stage; 
extension to other models is straightforward. 
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FIG. 1: QS order parameter vs. creation rate A in the CP; system sizes L = 10, 15 and 20 (upper 
to lower). Upper inset: moment ratio r^xx for system sizes 5, 10, 15, and 20, in order of increasing 
maximum value. Lower inset: entropy per site s for sizes 17, 19, 21, and 23 (left to right). 

I determined the QSD for the contact process on rings of up to 23 sites. The number 
of configurations scales as N c ~ (2 L — 2)/L + 1 (for L prime this formula is exact). The 
number of annihilation transitions is N a ~ 2 L ~ 1 (on average half the sites are occupied) and 
that of creation transitions is N cr ~ 3 ■ 2 L ~ 3 . (For L = 23, there are 364 723 configurations, 





FIG. 2: Contact process: values of A at crossings of r2ii (upper set), and at kurtosis minima (lower 
set) versus 1/L L5 . Inset: kurtosis q versus A for (lower to upper) L = 17, 19, 21, and 23. 

and ~ 7.3 x 10 6 transitions; the calculation takes about 16 hours on a 3 GHz processor.) 

The QS order parameter p (Fig. CD), follows the anticipated trend (i.e., p(A) is a sigmoidal 
function), but does not show any clear sign of the critical point; indeed, no such sign is 
expected for the small systems considered here. A precise estimate of the critical value A c 
can nevertheless be obtained through analysis of the moment ratios. Let rrij denote the 
j-th moment of the occupied site density, and r 211 = mijm^. The values A rj £, marking the 
crossing of r 2 ii(L) and r 2 n(L+l), approach the critical value systematically, as shown in Fig. 
EJ (A r _L is plotted, for convenience, versus L _L5 as this leads to an approximately linear plot.) 
Once preliminary estimates of X r> L (with uncertainty ~ 10~ 5 ) have been obtained, I perform 
high-resolution studies, with AA = 10~ 4 , in the vicinity of each crossing; precise estimates 
of the crossing values (uncertainty ~ 10~ 13 ) are then obtained applying Neville's algorithm 
[l9 ] to the data for r 2 n(L + 1) — r 2 n(L). Using the Bulirsch-Stoer (BST) extrapolation 
technique Q, 0], the data for sizes 8 to 23 furnish A c = 3.2961(15) and the critical moment 
ratio r2n iC = 1.1729(1). These values compare well with the best available estimates of 



A c = 3.297848(20) and r 2 ii jC = 1.1736(1) [13]; the associated errors are ~ 0.05%, remarkably 
small, in light of the system sizes used. 

The moment ratios r^m = m^/m\ also exhibit crossings that converge to A c . More 
surprisingly, the product m^imi exhibits crossings and appears to approach a well defined 
limit, 1.366(1), as L — > oo at the critical point; simulations (L = 1000 and 2000), yield 



1.374(2) for this quantity. The reduced fourth cumulant, or kurtosis, given by q(X, L) = 
K4/K2, where K2 = vn2 — m\ (the variance of the order parameter) and K4 = 1714 — 
4m 3 mi — 3m 2 , + 12m 2 m^ — 6m 4 , does not exhibit crossings but instead takes a pronounced 
minimum at a value \ qm (L) that converges to the critical value as A c — \ qrn (L) oc L -1,39 ^ 1 ). 
(This property has been verified in simulations using L = 1000, and 2000; departures from 



the minimum value are evident for |A| = 5 x 10—4 21].) The sharpness of the minimum, 
as gauged by q" = d 2 q/d\ 2 \\ qm , appears to increase rapidly with size: q" oc L 1 - 84 ^). (This 
is consistent with q" ~ L 2 ^ u± , as expected from FSS.) Since a negative kurtosis reflects a 
probability distribution that is broader at the maximum, and with shorter tails (compared 
to a Gaussian distribution with the same mean and variance), it is natural that q should be 
minimum at the critical point, where fluctuations are dominant. 

The statistical entropy per site, s = —L^ 1 J2jPj l n Pj; is plotted in Fig. 1. In the large-!/ 
limit, s should be zero for A < A c , since the QSD is concentrated on a set of configurations 
with vanishing density. As L —>■ 00, one expects ds/d\\\ c to diverge (as is the case for 
dp/dX), and to attain a maximum at some A > A c , approaching zero as A — > 00. The 
numerical data are consistent with these trends. 

Encouraged by the good results of the moment-ratio analysis, I examine three quantities 
expected, on the basis of FSS, to exhibit scaling at the critical point: the QS lifetime r, the 
order parameter, and the derivative r' = \dr<2ii/d\\. (The latter should diverge cx L 1 ^ u± .) 
Despite the small system sizes, these quantities indeed appear to follow power laws (see 
Fig. [3]). To obtain precise estimates of the associated exponents, I calculate the finite 
difference ratios Alnp/AlnL = [lnp(L) — p[L — l)]/[lnL — ln(L — 1)], (and similarly for 
t and r'). Linear regression of these ratios versus 1/L, using the seven largest sizes, yields 
P/v ± = 0.25193(3), z = 1.58054(2), and u ± = 1.092(1). These estimates differ by 0.06%, 
0.01%, and 0.5%, respectively, from the literature values of 0.25208, 1.5807, and 1.0968 [22]. 
Thus data on QS properties of systems with 23 sites or fewer yield estimates of critical 
exponents to within half a percent or better! Precise estimates are also found for r 3111 and 
q at the critical point: using BST extrapolation I find values of 1.5306(5) and -0.5015(5), 
compared with the simulation values of 1.526(3) and -0.505(3), respectively |23| . 

To test the robustness of this approach, I apply it to the pair contact process (PCP). 



In the PCP 



24j, each site is again either occupied or vacant, but all transitions involve 



a pair of particles occupying nearest-neighbor sites, called a pair in what follows. A pair 




FIG. 3: Main graph: finite differences Alnp/AlnL versus 1/L in the CP (upper set) and PCP 
(lower set); left inset: finite differences Alnri/AlnL (symbols as in main graph); right inset: order 
parameter in CP at A c versus system size, on log scales. 

annihilates itself at rate p, and with rate 1 — p creates a new particle at a randomly chosen 
site neighboring the pair, if this site is vacant. Any confi gura tion lacking a pair of nearest- 



neighbor occupied sites is absorbing. Simulation results [23j, |24j, |25[ place the PCP in the 
same universality class as the CP (namely, that of directed percolation). Unlike the CP, 
for which quite precise results have been derived via series expansions, there are no reliable 
predictions from series or other analytic methods. 

Using, as before, the parameter values associated with crossings of the moment ratio m^w, 
I obtain (for system sizes L = 8 to 23), the estimate p c = 0.07330(3), about 0.3% above 
the best available estimate of 0.077092(1) |26j]. The estimates @/v± = 0.2483(1) and v± = 
1.096(2), obtained via the same procedure as used for the CP, are also in good agreement 
with the accepted values. Analysis of the QS lifetime however, yields the unexpectedly 
large value z ~ 2.6. In fact, the finite-difference ratios Alnr/AlnL vary erratically with L, 
indicating that the result for z is unreliable. This may be associated with the large number of 
absorbing configurations in the PCP (growing exponentially with L), so that the extinction 
rate has not reached its asymptotic limiting behavior at the system sizes considered here. 



Extrapolation of the moment ratio at p c yields m 12 2 = 1.174(1), and reduced fourth cumulant 
q = —0.500(5), again in good accord with the expected values. (As in the case of the CP, the 
value of A at which q takes its minimum approaches A c with increasing system size.) Thus 
the QS properties of the PCP (for L < 23) permit one to assign the model to the directed 
percolation class, despite the lack of a clear result for the dynamic exponent z. 

It is natural to inquire whether the QS probability distribution exhibits any simplifying 
features. In an equilibrium lattice gas with interactions that do not extend beyond nearest 
neighbors, for example, the probability of a configuration depends only on the number 
of particles N and nearest-neighbor pairs P. The the CP, by contrast, I find that the QS 
probability of each configuration in a given (N, P) class is distinct (the probabilities typically 
vary over an order of magnitude or more, even far from the critical point). In a broad sense, 
this is because, unlike in equilibrium, not all annihilation events possess a complementary 
creation event. For similar reasons, it does not appear likely that the QSD could be obtained 
via the maximization of the statistical entropy, subject to some simple set of constraints. 

In summary, I show that analysis of exact (numerical) quasistationary properties on 
relatively small rings yields remarkably precise results for critical properties at an absorbing- 
state phase transition. Deriving the QS distribution involves rather modest programming 
and computational effort: the results reported here can be obtained in a few days on a fast 
microcomputer. Applied to the contact process, the analysis yields most critical properties 
with an error well below 0.1%. For the more complicated PCP, errors are generally < 
1%. Application to other absorbing-state phase transitions, including some belonging to 
other universality classes, is in progress. The method may also be useful in the study of 
metastable states, provides a valuable check on simulations, and may serve as the basis for 
phenomenological renormalization group approaches. 
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